suppressPackageStartupMessages({
  library(readr)
  library(tidyverse)
  library(reshape)
  library(popbio)
})

Qualitative description of population structure

Description of the demographic information

A life-cycle diagram for your species

DiagrammeR::grViz("
      digraph boxes_and_circles{

# Define nodes
node [shape = circle]

egg; a1; a2; a3; a4; a5; a6; a7; a8; a9; a10; a11; a12; a13

# Add edge statements
#Advaance classes
egg -> a1[label = m0]
a1 -> a2[label = m1]
a2 -> a3[label = m2]
a3 -> a4[label = m3]
a4 -> a5[label = m4]
a5 -> a6[label = m5]
a6 -> a7[label = m6]
a7 -> a8[label = m7]
a8 -> a9[label = m8]
a9 -> a10[label = m9]
a10 -> a11[label = m10]
a11 -> a12[label = m11]
a12 -> a13[label = m12]


# Fecundities
a1 -> egg[label = f1]
a2 -> egg[label = f2]
a3 -> egg[label = f3]
a4 -> egg[label = f4]
a5 -> egg[label = f5]
a6 -> egg[label = f6]
a7 -> egg[label = f7]
a8 -> egg[label = f8]
a9 -> egg[label = f9]
a10 -> egg[label = f10]
a11 -> egg[label = f11]
a12 -> egg[label = f12]
a13 -> egg[label = f13]
}
      ")

Will your model have a pre-breeding census or post-breeding census?

The matrix written out symbolically

Load the data

load("size.RData") #This is the entire data, loads it into a data.frame called size
# Define the functions we need

# Convert length to age using von bertalanffy model, solving for t
length2age <- function(length, l_inf, K, t_o){
  age <- (1/-K)*(log(1-(length/L_inf))) + t_o
  return(age)
}

# Convert age to length using von bertalanffy model
age2length <- function(age, l_inf, K, t_o){
  length <- l_inf*(1-exp(-K*(age-t_o)))
  return(length)
}

#Convert length to fecundity (number of eggs)
fecundity <- function(length, a, b){
  f <- 10^(a+(b*log10(length*10)))
  return(f)
}

Define demographic parameters

# Define Von Bert parameters
# Garbin & Castello, 2014
# L_inf <- mean(c(80, 80, 87.12, 94, 97.9, 97.3, 112.34, 89.38, 92.5))
# K <- mean(c(0.32, 0.6, 0.22, 0.38, 0.14, 0.25, 0.14, 0.38, 0.16))

# Or from Us and Stiurskjgfas, 1981
L_inf <- 102
K <- 0.55
t_0 <- -0.02

# Just to be clear, we use L-inf  and to from Ushisomething, and K from the review

# Define fecundity parameters
# From fishbase http://www.fishbase.se/Reproduction/FecundityList.php?ID=107&GenusName=Katsuwonus&SpeciesName=pelamis&fc=416&StockCode=121
fec_a <- -1.33354
fec_b <- 3.238
# Define mortality

m <- 0.63
z <- 1.39

Check population size by cohorts through time

size %>%
  mutate(Age = round(length2age(ClassFrq, L_inf, K, t_0))) %>%
  group_by(YearC, Age) %>%
  summarize(N = sum(as.numeric(Nr))) %>% 
  mutate(Age = as.factor(Age)) %>% 
  ungroup() %>% 
  ggplot(aes(x = YearC, y = log(N), color = Age)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  labs(x = "Year", y = "log N")

The matrix filled in with values

A <- matrix(0, 14, 14) #initial empty matrix with all 0

# Populate matrix with mortality
for (i in 2:14){
  A[i,i-1] <- exp(-1-z)
}

# Populate matrix with fecundity
ages <- seq(0:13)-1
lengths <- age2length(ages, L_inf, K, t_0)
A[1,] <- fecundity(lengths, fec_a, fec_b)

A[is.na(A)] <- 0
A[2,1] <- 0.0000001

colnames(A) <- ages
rownames(A) <- ages

knitr::kable(A, col.names = paste0("$a_",ages,"$"), row.names = F)
\(a_0\) \(a_1\) \(a_2\) \(a_3\) \(a_4\) \(a_5\) \(a_6\) \(a_7\) \(a_8\) \(a_9\) \(a_10\) \(a_11\) \(a_12\) \(a_13\)
114.4483696 1.657256e+07 7.026737e+07 1.294406e+08 1.758241e+08 2.072322e+08 2.27012e+08 2.38998e+08 2.461086e+08 2.502768e+08 2.527038e+08 2.541114e+08 2.549259e+08 255396708
0.0000001 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0.0000000 9.162970e-02 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0.0000000 0.000000e+00 9.162970e-02 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0.0000000 0.000000e+00 0.000000e+00 9.162970e-02 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 9.162970e-02 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 9.162970e-02 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 9.16297e-02 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 9.16297e-02 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 9.162970e-02 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0
0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 9.162970e-02 0.000000e+00 0.000000e+00 0.000000e+00 0
0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 9.162970e-02 0.000000e+00 0.000000e+00 0
0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 9.162970e-02 0.000000e+00 0
0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.00000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 9.162970e-02 0

Build a vector with 13 age classes

load("size_dist_1980.RData") #This is the n-vector for 1980, loads it to a vector called n

n <- n %>% {
  .$N}%>% 
  c(rep(0, times = 8))

n_0 <- sum(A[1, 2:14]*n, na.rm = T)

n <- c(n_0, n)

Project the population

project <- popbio::pop.projection(A, n, 30)
matplot(t(log(project$stage.vectors[-1,])), type = "l")